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Abstract 

Basic problems of the semiclassical microscopic modelling of strongly in- 
teracting systems are discussed within the framework of Quantum Molecular 
Dynamics (QMD). This model allows to study the influence of several types of 
nucleonic interactions on a large variety of observables and phenomena occur- 
ring in heavy ion collisions at relativistic energies. It is shown that the same 
predictions can be obtained with several - numerically completely different 
and independently written - programs as far as the same model parameters 
are employed and the same basic approximations are made. Many observ- 
ables are robust against variations of the details of the model assumptions 
used. Some of the physical results, however, depend also on rather technical 
parameters like the preparation of the initial configuration in phase space. 
This crucial problem is connected with the description of the ground state of 
single nuclei, which differs among the various approaches. An outlook to an 
improved molecular dynamics scheme for heavy ion collisions is given. 

25.75+r 
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I. INTRODUCTION 



One of the main interests of the study of relativistic heavy ion collisions is the investiga- 
tion of the properties of nuclear matter at extreme densities and excitation energies. . 
These investigations include the production of secondary particles, the properties of parti- 
cles in a (dense) nuclear medium, the compression and repulsion of dense nuclear matter, 
its equilibration during the reaction and its decay into fragments and single particles. On 
a macroscopic level the total energy of a dense nuclear system and its decomposition into 
thermal and compressional parts is related to the concept of the nuclear equation of state. 
Since a consistent derivation of the nuclear equation of state, e.g. the energy per nucleon as 
a function of density and temperature, is only possible in the low density limit (Bruckner 
theory) a reliable theoretical description is not at hand. On the other hand this quantity is 
of interest for many astrophysical questions [[14]] and therefore its knowledge is highly desir- 
able. Heavy ion reactions in combination with corresponding simulations using a variety of 
parametrizations of the equation of state are presently the only possible approach to study 
this quantity. 

Heavy ion collisions allow to search for a large number of observables which may be 
used as indicators of the properties of matter under extreme conditions. Frequently these 
observables are related to the quantitative description of collective effects like the bounce-off 
of cold spectator matter in the reaction plane M and the squeeze-out of hot and compressed 



participant matter perpendicular to the reaction plane [ID as well as to the production of 
secondary particles [PH 131 . 

Experiments performed at LBL in the early 80's (Streamer chamber, Plastic ball) yield 
first 4tt information of the final momentum distributions in heavy ion reactions [|l^,|nj . New 
experimental 4n setups at LBL, Ganil, GSI and Brookhaven enable precise measurements 
on the emission of primary and secondary particles and therefore provide a stimulating 
challenge to the theoretical description of heavy ion collisions 

Lots of comparisons have been made between experimental data and microscopic and 
macroscopic transport-theoretical calculations. Besides other microscopic models like VUU 
T7J, BUU 0,0, Landau- Vlasov AMD ^ or FMD the Quantum Molecular Dy- 



namics approach (QMD) is a frequently used model [23]. However, from recent comparisons 



of experimental results with QMD using different numerical realizations conflicting results 
have been reported pl| . We will demonstrate that these discrepancies are on the one hand 
due to the variation of physical parameters (like ground state densities, interaction ranges) 
whose precise values are not known. On the other hand they are a consequence of the im- 
possibility to build a ground state nucleus with all its detailed structure in a semiclassical 
molecular approach. 

This paper is organized as follows: First the basic principles of a microscopic modelling 
of heavy ion reactions are briefly reported. The assumptions entering in the different QMD 
realizations are described in detail. The origin of differences is critically examined in this 
context. We demonstrate that most of the discrepancies can be attributed to different 
descriptions of the initial nuclei, which limits the applicability of some versions. Finally an 
outlook to a new molecular dynamics scheme for heavy ion collisions simulations is given. 
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II. MICROSCOPIC MODELLING OF HEAVY ION REACTIONS 



Presently the microscopic models can be subdivided into two classes: Those which follow 
the time evolution of the one-body phase space distribution and those which are based on 
n-body molecular dynamics or cascade schemes. 



A. VUU-type models 

The microscopic transport models for the one-body Wigner phase space density distribu- 
tion obtained different names although they solve the same equation. They differ in the tech- 
nical realization, i.e. the computer program, and are known as Vlasov-Uehling-Uhlenbeck 
(VUU) model [0,0 (or BUU [0,0, LV J20| ). They solve the following transport equation 
for the one-body Wigner density f(f,p,t) in the limit h — > 0: 

df,^-, T-j TTX7f 47T 3 (ftc) 4 f dV x d s p' 2 j3 da 

- + v .V r f-V r U.V p f = -—^j—— dp2m 

x[fh{i-f[){i-f' 2 )-f[m-f){i-h)\ 

X S A (p + p 2 -p[ -p' 2 ). (1) 

The l.h.s. of this equation is the total differential of / with respect to the time assuming 
a momentum independent potential U. This potential is calculated selfconsistently and 
corresponds to the real part of the Bruckner G-matrix. Usually a Skyrme-parametrization 

of the real part of the G-matrix is employed, where p is the nuclear density which is frequently 
measured in units of the saturation density po of cold nuclear matter. 

The r.h.s. of Eq. (|l|) contains a Boltzmann collision integral, which is identified with 
the imaginary part of the G-matrix. This part describes the influence of binary hard-core 
collisions, where the term with ff'2 describes the loss of particles (in a phase space region) 
and the term with f[f 2 the gain term due to collisions feeding the considered phase space 
region. It is supplemented with the Nordheim-Uehling-Uhlenbeck modifications in order 
to obey the Pauli-principle in the final state of the collisions [T2BJ. The (^-functions assure 
the conservation of the four- momentum. The cross section a is normally adjusted to the 
free nucleon-nucleon scattering. The differences from cross sections calculated from the 



imaginary part of the Bruckner G-matrix are minor |5T| and influence little the observables 
of a heavy ion collision. For a derivation of this equation see [44 , 45] , 

The equation is solved by use of the testparticle method. Here the continuous one-body 
distribution function / at t — is represented by an ensemble of n - (A p + A t ) pointlike parti- 
cles. This is often viewed as an ensemble of n parallel events with A p + A t physical particles 
each, where A p and A t denote the number of nucleons in projectile and target, respectively. 
The l.h.s. of Eq. (]!]) can be regarded as the transport equation (Vlasov-equation) for a dis- 
tribution of classical particles whose time evolution is governed by Hamilton's equations of 
motion. 
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(3) 



The testparticles move due to their own, selfconsistently generated mean-field. The r.h.s. 
is taken into account by additional stochastic scattering similar to the collisions in cascade 
models 0,0. 

More explicitly the test particle method corresponds to the replacement of the expecta- 
tion value of a single particle observable 



(0(t)) = J f(r,p,t)0(r,p)d 3 rd 3 p (4) 
by a Monte Carlo integration 

, n(A T +A P ) 

m))= — A £ om)Mt)) (5) 

n{A T + A P ) . =1 

where the fi(t) and Pi(t) are points in phase space which are distributed according to 
f(p,r,t), i.e., 

j n(A T +A P ) 

f(p,r,t)= lim £ <5(f-rl(t))<5(p-pl(t)) (6) 

n-^oo n (A T + Ap) g 

It is evident that a large number n is necessary to avoid numerical noise. Predictions 
beyond the one-body level are not feasible although several attempts have been made to 
relate the (unphysical) numerical noise to physical fluctuations. In practice the number n 
lies in the range between 15 and 500 and one employs a grid to obtain a smooth phase space 
density distribution. 

The numerical realization can be achieved in various ways. VUU uses a phase space 
sphere around each particle in order to determine / and a coordinate space sphere to deter- 
mine p and thus U(p). This corresponds to a Lagrangian method. On the contrary, BUU 
uses a fixed grid corresponding to an Eulerian method in hydrodynamics. In both models 
collisions are treated in a parallel event method, only testparticles of the same events, i.e. 
the A p + At test particles with the same index n, can collide. The Landau- Vlasov model 
determines / by the overlap of several Gaussians. The collisions are performed in a crossed- 
event (or full ensemble) method where all n(A p + A t ) may collide with each other particle 
with a scaled cross section. 

For a solution of Eq. (|I|) proper boundary conditions have to be specified. In the case of 
heavy ion reactions, the test particles are distributed according to the density- and (Fermi-) 
momentum distribution of ground state nuclei. The latter are then boosted onto every other 
with the proper relative momentum. Initially the test particles are randomly distributed in 
a coordinate space sphere of the radius R = 1.12A 1 / 3 fm (where A is the atomic number 
of the nucleus) and in a momentum space sphere of the radius of the corresponding Fermi 
momentum. 

One should keep in mind that the forces acting on the testparticles are calculated from 
the entire distribution including testparticles from all events, hence the n parallel events 
are not independent and event-by-event correlations cannot be analyzed within this one- 
body transport models. In the limit n — > oo the distribution of these propagated test 
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particles at the time t represents the one-body distribution function at this time. Any one- 
body observable can be calculated by averaging the values weighted with the distribution 
function according to Eq. (|5|). Hence, VUU type models succeeded in the description of one- 
body observables like collective flow, stopping and particle spectra, but, fluctuations and 
correlations, such as the formation of fragments or the description of two-particle correlations 
in relativistic heavy ion collisions, are beyond the scope of a transport model based on 
a one-body distribution function Any fluctuation of the observables seen in the 

Monte Carlo simulation of the one-body distribution function is due to numerical noise and 
disappears in the limit of a infinite number of test particles. 

B. The Quantum Molecular Dynamics approach 

An approach which goes beyond a one-body description as explained above, is the Quan- 
tum Molecular Dynamics (QMD) model EI,E3,B^,H] ■ The QMD model is a ra-body theory 
which simulates heavy ion reactions at intermediate energies on an event by event basis. 
Taking into account all fluctuations and correlations has basically two advantages: i) many- 
body processes, in particular the formation of complex fragments are explicitly treated and 
ii) the model allows for an event-by-event analysis of heavy ion reactions similar to the 
methods which are used for the analysis of exclusive high acceptance data. 

The major aspects of the formulation of QMD will now be discussed briefly. For a more 
detailed description we refer to ref. ||. The particular realizations of this model will be 
discussed later. 

1. Formal derivation of the transport equation 

In QMD each nucleon is represented by a coherent state of the form (we set H,c — 1) 
which are characterized by 6 time-dependent parameters, r*j and Pi, respectively. 

/ O \3/4 

( j ) .(£, t )=l — \ e -{^i-m) 2 /L e i^iPiit)_ ( 7 ) 

The parameter L, which is related to the extension of the wave packet in phase space, is 
fixed. The total n-body wave function is assumed to be the direct product of coherent states 

(0) 

$ = Y[<l>i(x i ,r i ,p i ,t) (8) 

i 

Note that we do not use a Slater determinant (with (A p + A t )\ summation terms) and 
thus neglect antisymmetrization. First successful attempts to simulate heavy ion reactions 
with antisymmetrized states have been performed for small systems |2l|j2^|. A consistent 
derivation of the QMD equations of motion for the wave function under the influence of 
both, the real and the imaginary part of the G-matrix, however, has not yet been achieved. 
Therefore we will add the imaginary part as a cross section and treat them as in the cascade 
approach. How to incorporate cross sections into a antisymmetrized molecular dynamics is 
not yet known. This limits its applicability to very low beam energies. 
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The initial values of the parameters are chosen in a way that the ensemble of At + Ap 
nucleons gives a proper density distribution as well as a proper momentum distribution of 
the projectile and target nuclei. 

The equations of motion of the many-body system is calculated by means of a generalized 



variational principle: we start out from the action [37 



S = J C[$,$*]dt 

ti 



(9) 



with the Lagrange functional £ 









ih--H 


dt 



(10) 



where the total time derivative includes the derivation with respect to the parameters. The 
Hamiltonian H contains a kinetic term and mutual interactions Vij, which can be interpreted 
as the real part of the Bruckner G-matrix supplemented by the Coulomb interaction. We 
will lateron describe the components of H in detail. The time evolution of the parameters 
is obtained by the requirement that the action is stationary under the allowed variation of 
the wave function. This yields an Euler-Lagrange equation for each parameter. 

If the true solution of the Schrodinger equation is contained in the restricted set of 
wave functions (f>i(xi,t) (with parameters fi,Pi) this variation of the action gives the exact 
solution of the Schrodinger equation. If the parameter space is too restricted we obtain 
that wave function in the restricted parameter space which comes closest to the solution of 
the Schrodinger equation. Note that the set of wave functions which can be covered with 
special parametrizations is not necessarily a subspace of Hilbert-space, thus the superposition 
principle does not hold. 

For the coherent states and a Hamiltonian of the form H = ^ i T i + | YUj Vij kinetic 
energy, = potential energy) the Lagrangian and the variation can easily be calculated 
and we obtain: 



2Lm 



(11) 



Variation yields: 



fit 



(12) 



a = -v«i;(fy) = -v*<#) (is) 

with (Vij) = J d 3 x\d 3 X2(f)*(j)*jV(xi,X2)4 > i ( Pj- These are the time evolution equations which 
are solved numerically. Thus the variational principle reduces the time evolution of the 
n-body Schrodinger equation to the time evolution equations of 6 • (Ap + At) parameters to 
which a physical meaning can be attributed. The equations of motion for the parameters Pi 
and r*j read 
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(14) 



and show the same structure as the classical Hamilton equations, Eq. (|3]). The numerical 
solution can be treated in a similar manner as it is done in classical molecular dynamics 
|33H36| . Trial wave functions other than the gaussians in Eq. (0) , yield more complex equa- 
tions of motion for other parameters and hence the analogy to classical molecular dynamics 
is lost. If (H) has no explicit time dependence, QMD conserves energy and momentum by 
construction. 



2. Description of the Hamiltonian 

The nuclear dynamics of the QMD can also be translated into a semiclassical scheme. 
The Wigner distribution function of the nucleon i can be easily derived from the test wave 
functions (note that antisymmetrization is neglected). 
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fi(r,P,t) = ^e-^-^ie-^^ 2 ^ (15) 

and the total Wigner density is the sum of those of all nucleons. Hence the expectation 
value of the total Hamiltonian reads 

(H) = (T) + (V) 

2 

= E^ + EE/ W^V^f^r'J'^drdr'dpdf . (16) 

i 1 i j>i 

The baryon-potential consists of the real part of the G-Matrix which is supplemented by the 
Coulomb interaction between the charged particles. The former can be further subdivided 
in a part containing the contact Skyrme-type interaction only, a contribution due to a finite 
range Yukawa-potential, and a momentum dependent part. V* J ' = G^ + Vy uk + V^ oul + V^ di 
consists of 



^Skyrme + ^Yuk + 'mdi + ^Coul 



ttffc - x 3 ) + Mixi - x.v- 1 ^) + t3 ^_K_MM + (17) 

I *^ I / r 1 

ZZe 2 

t 4 ln 2 (l + h{pi - Pj) 2 )S{xi - xj) + ^ J 



Zi, Zj are the charges of the baryons % and j. The real part of the Bruckner G-matrix is 
density dependent, which is reflected in the expression for G l K The expectation value of G 
for the nucleon % is a function of the interaction density p\ nt . It is indeed this quantity which 
relates the number density to the energy content of nuclear matter. 

PL(fd= (^ge-W-^ (18) 
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Note that the interaction density has twice the width of the single particle density. Moreover, 
the particles do not interact with themselves. This is different compared to VUU-type models 
because in QMD explicit N-N interactions are treated, hence the force acting on a particle 
at the position r depends on the exact positions of all other particles, whereas the density 
employed in the one-body theories (eq. (|6])) depends on the average number of nucleons in 
the vicinity of the test particle only. 

It should be noted that the width L of the distribution function determines the interaction 
range of the particle and influences the density distribution of finite systems. Therefore its 
value has to be adopted to reasonable interaction ranges of the strong interaction. 

The momentum dependence V^ di of the N-N interaction, which may optionally be used 
in QMD, is fitted to experimental data |38|]3]| on the real part of the nucleon optical potential 



BUggi, which yields 

U mdl = 5 • In 2 (e ■ (Ap) 2 + l) • (19) 

These measurements have been superseded recently by new data [42] and thus a new 
parametrization has been advanced ^3] . 



The potential part of the equation of state (we will discuss this concept in the next 
subsection in more detail) resulting from the convolution of the distribution functions fi and 
fj with the interactions V^ yime + (local interactions including momentum dependence) 
then reads: 

U = a- + (3 ■ (^j 7 + 5 • In 2 (e ■ (Ap) 2 + l) ■ (^j (20) 

Here it should be noted that due to the definition of p int (eq. [TJ) no mean-field potentials 
(as e.g. eq. (0) for VUU) show up in the calculation of the equations of motion (eq. [14] ) 



of QMD but a sum of two (and three) body interactions (see eq. |TTD . Hence energy and 
momentum are - in contrast to single VUU 'events'- strictly conserved in each event. 

The Coulomb interaction cannot be treated for infinite matter, since this leads to diverg- 
ing terms. In the first versions of QMD no explicit treatment of the isospin is performed 
and the charges are replaced by effective charges, i.e. all nucleons had been attributed the 
effective charge Z = [Z vro ^ + Z targ ) / (A pro j + A targ ). IQMD (we will later come to that) 
and other more recent versions use the real baryon charges. 

The parameters /i and t\...t$ are adjusted to fit the real part of the G-matrix and to 
describe the properties of finite nuclei. 



3. The relation to the nuclear equation of state 

One strong motivation for the numerical simulation of heavy ion reactions is the possi- 
bility to investigate effects of the underlying nuclear equation of state on the dynamics and 
final states of these collisions. QMD is a model for non-equilibrium dynamics with mutual 
interactions among the constituents and therefore does not contain any parametrization of 
the nuclear equation of state in terms of an explicit relation between number density, tem- 
perature and the energy density. In equilibrium and in the thermodynamic limit (n — ► oo), 
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however, such a functional relation can be deduced from the nucleon-nucleon potentials and 
the cross-sections employed in the model. 

For the description of the energy per nucleon as a function of density (assuming T = 0) 
usually Skyrme type parametrizations (see eq. (||)) are used. This ansatz is phenomenolog- 
ical and can be derived for the case 7 = 2 from the assumption that the particles interact 
with each other with two- and threebody contact forces. It is generalized to effective higher 
order contact terms by setting 7 > 1 to be a real number. This generalized ansatz uses three 
parameters a, /3, 7; two of them are fixed by the constraint that the total energy should have 
a minimum at the saturation density p = po with a value of E/A = — 16MeV which cor- 
responds to the the volume energy in the Bethe-Weizsacker mass formula. Together with 
the condition that a free particle has no binding energy (which is automatically fulfilled 
within this ansatz) there remains one degree of freedom. The third parameter is fixed by 
the nuclear compressibility, which is the second derivative of the energy at the minimum 
with respect to the density: 

-VM) (21) 



dp 2 \A 

Two different equations of state are commonly used: A hard equation of state (H) with 
a compressibility of k =380 MeV and a soft equation of state (S) with a compressibility of 
k =200 MeV [T7y27[. 

To derive an equation of state from the interactions used in eq. (fl7|) we have to convolute 
the potentials with the distribution functions assuming an infinite homogeneous distribution. 
In this limit the Vskyrme and Vyuk become functions of the constant density only. The 
interaction density of eq. ( pT8|) as used in eq. (|20| ) can be replaced by the position independent 
nuclear matter density. The integration over the relative momenta of infinite nuclear matter 
Fermi distributions finally turns into a density dependence of the momentum dependent 
interaction. This allows us to obtain the compressional part of the nuclear equation of state, 
which depends on the density only. The parameters of the interactions in eq. flT7| ) can 
therefore be chosen that way that a hard or soft eos is obtained for the infinite matter case. 
It should again be noted that the parameters of the potentials allow a relation to the nuclear 
equation of state (eos) but that the microscopic description works as well for systems far off 
from equilibrium where no eos can be defined. 

The interaction range parameter L influences the interaction density (eq. |I3|) for finite 
systems. For (homogeneous) infinite nuclear matter the density (and thus the potential 
energy) do not depend anymore on the extension of the gaussian wavepackets. Thus, the 
equation of state of infinite nuclear matter is independent of L. In finite matter E/A also 
depends on L. Thus even two parametrizations which yield the same eos may produce 
different results for the reaction of two heavy ions. Therefore we have to adjust L to have 
reasonable surface properties. In order to allow a physical interpretation L should be in the 
order of the size one expects for the range of the nuclear interaction. There exists a range 
of values for L, which allows to fix these properties. Larger values of L increase the effective 
range of the interaction and thus lead to some smearing of fluctuations, which are stronger 
for more located wavepackets (small values of L). 

Hence, the nuclear equation of state can only be defined as the bulk properties in the 
limit of an infinite system: The concept of the nuclear equation of state as discussed here 
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does only make sense for large macroscopic systems in (at least local) equilibrium, while the 
ansatz with mutual interactions has no restrictions with respect to the size and is therefore 
also applicable for finite systems far off equilibrium. The time evolution of the non-thermal 
system of two reacting heavy ions is completely determined by the two-body potentials and 
the scattering cross sections, respectively. 

In QMD the parameters t\...t$ are uniquely related to the corresponding values of 
a, [3, 7, 5 and e which serve as input. The standard values of these parameters can be 
found in table |. 



4- Inclusion of collisions 

As stated above the imaginary part of the G-matrix acts like a collision term. In the 
QMD simulation we restrict ourselves to binary collisions (two-body level). The collisions 
are performed in a point-particle sense in a similar way as in VUU or cascade: Two particles 
collide if their minimum distance d, i.e. the minimum relative distance of the centroids of 
the Gaussians during their motion, in their CM frame fulfills the requirement: 



d<d = J—, otot = a[\/s, type). (22) 

where the cross section is assumed to be the free cross section of the regarded collision type 
(N-N, N- A, ...). 

Beside the parameters describing the N-N potential, the cross sections constitute an- 
other major part of the model. In principle, both sections of parameters are connected 
and can be deduced from Bruckner theory. QMD-calculations using consistently derived 
cross-sections and potentials from the local phase space distributions have been discussed 
e.g. in [72| . Such simulations are time-consuming since the cross-sections and potentials do 



explicitly depend on the local phase space population. 

Within the framework of using free cross section one may parametrize the cross section 
of the processes to fit to the experimental data if available. For unknown cross sections 
isospin symmetry and time reversibility is assumed. 

Alternatively, cross-sections may be obtained from theoretical considerations. For one 
particular QMD- version the one boson exchange model has been employed for this purpose. 
This has the advantage to have a first handle for the description of cross sections in the 
nuclear medium. 

If two particles scatter, the direction of the final momenta will be distributed randomly 
in such a way that the distribution of many identical collisions corresponds to the measured 
cross section. For elastic scattering the distribution is taken from [j47| : 



^~exp(A( S )-t) , (23) 

where t is — q 2 , the squared momentum transfer (which also includes the information on the 
polar angle) and y/s is the cm. energy in GeV. 

It should be noted that the presented treatment of the collisions may cause problems 
with causality since the particles can interact immediately at a distance. The collision 
information is given to both particles at the same time when they are at closest position. 
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It should also be noted that the time order of the collisions is determined in a common 
system of all particles. The evolution of the system is propagated with one common clock. 
As it has been already pointed out by Kodama et al. ||8| the time ordering is not unique. 
Thus the choice of the common referential system may influence the observables. Normally 
a system is chosen where the relative velocities with respect to that system are as small as 
possible. Thus BQMD used the nucleus-nucleus CM system as referential system and VUU 
and IQMD use the nucleon-nucleon CM system. The choice of the Lab system as referential 
system would e.g. cause for the system Au(lAGeV)+Au b=3fm, hard eos, an enhancement 
of the flow (in IQMD pf r rises from 98 ± 3 MeV/c to 110 ± 3 MeV/c in the Lab system ) 
and a reduction of the pion number (in IQMD N n falls from 64 ± 1 to 60 ± 1 in the Lab 
system) . 

Also the choice of the minimum distance point as collision point can be motivated within 
this respect. An earlier collision (e.g. at the point when the distance is sufficient to fulfill 
the distance condition) could cause stronger acausalities. It will also reduce the mean free 



path and thus enlarge stopping and flow [49 



5. Pauli blocking due to Fermi statistics 

The cross section is reduced to an effective cross section by the Pauli-blo eking. For each 
collision the phase space densities in the final states are checked in order to assure that the 
final distribution in phase space is in agreement with the Pauli principle (/ < 1). Phase 
space in QMD is not discretized into elementary cells as in one-body models like VUU, 
in order to obtain smooth distribution functions the following procedure is applied: The 
phase space density f[ at the final states 1' and 2' is measured and interpreted as a blocking 
probability. Thus, the collision is only allowed with a probability of (1 — f[)(l — fy). If the 
collision is not allowed the particles remain at their original momenta. 

The Pauli blockers of VUU and QMD show efficiencies of about 94-96 %, i.e. a single 
ground state nucleus with Fermi momentum would show a blocking rate of this amount. In 
order to reduce the noise of spurious collisions in ground state nuclei additional conditions 
allow a nucleon only to collide with a nucleon of the other nucleus or with a nucleon that has 
already undergone a collision. Nevertheless the problem of Pauli blocking causes a limitation 
of the calculated system to have not less incident energy than about the Fermi energy. 



C. Numerical structure 

The QMD model consists of three major parts, namely i) the initialisation of projectile 
and target, ii) the propagation of nucleons, resonances and newly produced particles due 
to their mutual potential interactions, and iii) the hard collisions according to the energy 
dependent cross section for the various channels together with the Pauli-blocking. 

For the propagation the description of the potential (or to be more exact of the real 
part of the Bruckner G-matrix) is of crucial importance. 

The solution of the transport equations for the iV-body distribution function is done in 
the following way: 
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1. Projectile and target are initialized. For each of these nuclei the nucleons initialized 
according to a distribution f(r,p,t = 0). This distribution is essentially constrained 
by the requirement to reproduce the ground state properties of the two nuclei, i.e. 
radii, binding energies. 

2. The particles are propagated using Hamilton's equations of motion (|TJ]) with a given 
Hamiltonian (H). 

3. Two particles close in coordinate space may perform a collision. The particles change 
their momenta respecting the Pauli principle. 

The input into the program may be subdivided into three classes of parameters 

Reaction parameters: projectile and target masses (and charges), bombarding energy, 
impact parameter. They define the whole kinematics of a single event. 

Physics Parameters: interaction range, potential parameters, in medium cross sections 
and decay widths, etc. They correspond to a detailed description of interactions and 
may be changed within a reasonable range. Finally their deduction is a particular goal 
of the comparison between calculation and experiment. 

Technical parameters: time step size, initial distance, cutoff parameters, maximum col- 
lision distance, etc. They are used to perform effective calculations on a computer. 
The observables should not depend on them. 

If all these parameters are fixed the calculation of a single event can be performed in the 
following way: 

• initialize projectile and target nuclei in their "ground states" as mentioned above, 

• propagate the constituents of the system according to their mutual potential and hard 
scattering interactions, this includes 

— calculation of interaction densities, forces and the Hamiltonian 

— propagation of all particles according to Hamilton's equation of motion 

— perform all collisions within this time step. Decide for each collision whether its 
final state is Pauli-blocked. If this is the case: keep momenta of collision partners 
unchanged, otherwise change momenta according to the angular distribution of 
this particular channel. 

• output of information (coordinates, momenta, scattering partners, . . . ) about the 
intermediate reaction stages and output of the final phase-space configuration (which 
would correspond to the freeze-out configuration in a thermal picture). 

This procedure is repeated until sufficient statistics, i.e. a large number of independent 
events, is obtained. 

This principal structure is common for all QMD realizations, which differ, however, in 
details and the initialisation of projectile and target. In the following we will study the 
influence of these differences on observable and nonobservable quantities. 
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III. DESCRIPTION OF PARTICULAR QMD MODEL REALIZATIONS 



The original QMD [0,^] program was developed further to include momentum depen- 



dent interactions pl] , |32| 



A. BQMD 



The original QMD has been rewritten by Bohnet et al. |5T[] for the purpose of studying 
low energy fragmentation data. This program has been dubbed BQMD since it was designed 
for describing the proper binding of a nucleus in order to describe fragmentation processes 
52-^5 . An improvement on the stability against artificial particle evaporation has been 



achieved in BQMD by a procedure explained below which causes fluctuations of the energy 
around the mean value by 2 MeV/nucleon. 



1. Initialisation in BQMD 

In BQMD the nucleons are distributed within a sphere with a Wood-Saxon-type den- 
sity profile. (The original QMD used a sphere for the distribution of the centroids of the 
Gaussians.) The maximum Fermi-momentum is limited by the local binding energy of the 
nucleon in order to keep all particles bound. By this procedure, however, the mean kinetic 
energy of the particles is lowered to about 10 to 12 MeV/nucleon. The ground state cen- 
tral density is assumed to be po = 0.15/m~ 3 . The Gaussian width for the interactions are 
chosen to be L = 4.33 fm 2 . The binding energy as given in the Weizsaecker mass formula 
is reproduced from Lithium up to the heaviest nuclei ||. As already seen in fig. 10 of || 
this particular version suffers from fluctuations of the rms radius. The consequences will be 
discussed later. 



2. Potentials in BQMD 

The range and the strength of the Yukawa potential in BQMD has been chosen to 
describe the surface of the nucleus best. In order to keep the nuclear equation of state and 
the binding energy independent of the Yukawa interactions and to keep the binding energy 
at its experimental value, the coupling constant t\ of the Skyrme-type two body interaction, 
is modified according to || 

t \(i)p(r7 ) = t lP (r7 )-Y< U Z Uk - ( 24 ) 

j 

Note that the Skyrme and Yukawa coupling constants are different for each particle here. 
With this procedure the validity of Newtons theorem actio = reactio can be assured on the 
ensemble average only, which also leads to violation of energy conservation in single events. 
The energy fluctuates about 2 MeV/nucleon around the mean value H. The range of the 
Yukawa-potential is chosen as 1.5 fm. 
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3. Collision term in BQMD 



BQMD has in common with the original QMD that it uses nucleons and deltas only. 
The employed cross sections have been parametrized by Cugnon fl4"? |. All nucleons interact 
with the same average cross section without distinction in isospin. The elastic cross section 
is given by a constant value of 55 mb for collisions with yfs < 1.8993 GeV and for higher 
energies by the parametrization: 

35 

<7ei(mb) = ; c + 20 ( 25 ) 

v ; 1 + 100 • (-y/s/GeV - 1.8993) v ' 

The inelastic cross section NN — > NA is zero for yfs < 2.015 GeV and for higher values by 
the parametrisation 

20r 2 

a m (mb) = — x = ^s/GeV - 2.015 (26) 

u. id oc 

da i 

The angular distribution of the collisions is described by ~ exp(A(s) t) with 
() _ (3.65 (VS/GeV - 1.8766)) 6 

1 ; l + (3.65(v^/GeV- 1.8766)) 6 ' 1 } 



B. IQMD 



The Isospin-QMD (IQMD) [0,^] treats the different charge states of nucleons, deltas 
and pions explicitly, as inherited from the VUU model. IQMD has been used for the analysis 
of collective flow effects of nucleons [p^,|53]-|63| and pions pl|-|B"2"||. Comparisons to exper- 
imental data with this model have been presented in [|66| |68|| . As it has been developed 



from the VUU-model, its coding is therefore independent of the original QMD. The isospin 
degrees of freedom enter into the cross sections (here cross sections of VUU |I7| similar to 
the parametrizations of VerWest and Arndt |59| have been taken, see also ref. |62 



as well 

as in the Coulomb interactions. The elastic and inelastic cross sections for proton-proton 
and proton-neutron collisions used in IQMD are shown in figure |l[ The cross section for 
neutron-neutron collisions are assumed to be equal to the proton-proton cross sections. 



1. Potentials used in IQMD 

The IQMD-model offers rather stable density distributions and good energy conservation, 
however for the price of nucleon evaporation and and improper binding energies {Eu n d ~ 4—5 
MeV/nucleon for heavy nuclei instead of 8 MeV/nucleon). 

In addition to the use of the explicit charge states of all baryons and mesons a symmetry 
potential between protons and neutrons corresponding to the Bethe-Weizsacker mass formula 
has been included 

VSm = tB-TxTySft - rj) t 6 = 100 MeV (28) 
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where T 3i and T 3 j denote the isospin projections of particles i and j. Other baryonic po- 
tentials like V S k yrme and V^ di are defined isospin-independent like in all other flavors. The 
Yukawa potential in IQMD is very short ranged (p = 0.4fm in contrast to p = 1.5 fm 
in BQMD) and weak. The modification of the a term of the static potential is done in 
an particle independent way. As in BQMD this corresponds to the interpretation that an 
additional term in the Skyrme ansatz which is proportional to (Vp) 2 can be expanded in 
first order to a term linear in density (which reduces a effectively) plus Yukawa potentials. 
Additional attractive Yukawa forces hence modify the EOS (and therefore the a term has 
to be modified to obtain the same EOS). Yukawa forces stabilize the nuclei because of the 
increase of the interaction range as compared to a 5-like Skyrme-potentials. Thus nucleons 
notice earlier that they will arrive at the surface and are more effectively decelerated as 
without this potential. In addition the fluctuations are reduced. 



2. Pions in IQMD 

Free pions are moving under the influence of the Coulomb interactions. Pions may be 
produced by the decay of a A-resonance and may be reabsorbed by a nucleon forming a delta 
again. IQMD and HQMD, which will be described in the next section, differ concerning the 
pion production in the production cross sections (HQMD uses cross sections based on the 
one boson exchange model), the included resonances (HQMD contains additionally iV* and 
NN — ► AA collisions) and the angular distribution of inelastic collisions (HQMD has more 
realistic non-isotropic distributions obtained from OBE calculations which are not present 
in original IQMD). Recent updates of IQMD calculating the pion production (e.g. [ |6"T| , |62"| ) 



also use the inelastic angular distributions of HQMD. The effect of this modification on 
nucleonic observables is quite small. 



3. Initialisation in IQMD 

The most important difference to BQMD is the initialisation. In IQMD the centroids of 
the Gaussians in a nucleus are randomly distributed in a phase space sphere (r < R and 
V ^ Pf) with R = A 1 / 3 ■ 1.12 fm corresponding to a ground state density of po — 0.17/mr 3 . 
The Fermi momentum pp depends on the ground state density. For po = 0.17/m -3 it has 
a value of about pp ~ 268 MeV/c. While, as said, in BQMD the maximum momentum 
is determined by the local binding energy (which causes an effective reduction of the total 
Fermi energy to about 10 - 12 MeV), in IQMD the momenta are uniformly distributed within 
a momentum sphere p < pp erm % ~ 268MeV/c without further local constraints. Therefore 
it may happen that nucleons close to the surface, where the local potential energy is low, 
are unbound initially. This possibility is not given in BQMD or HQMD. It gives, however, 
a reduced binding energy per nucleon as compared to the Weizsacker mass formula. Hence 
the initialized nuclei are less stable against spurious particle evaporation as compared to 
BQMD. On the other hand this ansatz makes available the full Fermi-energy calculated 
from the Skyrme ansatz. The full Fermi pressure yields (as compared to BQMD) a stronger 
stability of the density profile against vibration modes. Finally it should be noted that 



15 



IQMD performs a Lorentz contraction of the nucleus coordinate distribution which is not 
present in BQMD and which becomes important for higher energies E / nucleon > 1 GeV. 



4- Interaction range 

As it has already been stated, the Gaussian width can be regarded as a description of the 
interaction range of a particle. Its influence disappears for infinite nuclear matter whereas 
for finite systems it may play a non negligible role. 

In IQMD the Gaussian width can be used as an optional input parameter. The default 
version of uses a system dependent Gaussian width while BQMD uses L = 4.33 fm? inde- 
pendent of the system size. The system dependence of L in IQMD has been introduced 
in order to obtain maximum stability of the nucleonic density profiles. As an example for 
Au+Au a value of L = 8.66 fm 2 is choosen, for Ca+Ca and lighter nuclei L = 4.33. 



C. HQMD 

HQMD is an upgrade of QMD which combines optional features of BQMD and IQMD. 
It does not remedy the shortcomings of BQMD and IQMD, but allows to study the influence 
of the different modules on physics results. In addition, higher resonances ( the AT* (1440)), 
free pions and the proper isospin coupling have been incorporated by Huber et al. |)6| . The 
isospin degrees of freedom play an important role especially for the particle production. The 
employed inelastic cross sections NN — > iViV* , NA and A A have been calculated within an 
one-boson exchange model (OBE). Also the angular distribution of the inelastic reactions 
was calculated and parametrized in the following way: 

exp(6(s) -costf) , (29) 

a(s) and b(s) are functions of \fs and vary in their definition for different intervals of yfs (see 
table 0). is the polar angle. It should be noted that VUU, BUU and IQMD (in its older 
version) assumed isotropic scattering for the inelastic channels which causes differences in 
the flow at higher energies. For elastic collisions a new parametrisation |57J has been used 



which can be taken from table III 



In addition the numerical propagation routines have been changed to a higher accuracy. 
A 4th order Runge-Kutta propagation scheme allows an energy conservation of about 1 per 
mille. The Yukawa interaction has been suppressed. 

This upgrade of QMD which has been dedicated to the question concerning the meson 
production. It is quoted as HQMD because it contains higher resonances. 

HQMD offers the possibility to choose between the two initialisation modes of BQMD 
and IQMD. Moreover one can choose between the different parametrizations of the cross 
section as described above: the cross section parametrization used in IQMD and that used 
in BQMD. It was checked that it reproduces the results obtained with BQMD and with 
IQMD if the corresponding subroutines are used. Therefore it may be directly used to 
analyse the effects of the different ingredients in the QMD flavours. 
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D. Other flavours 



There exist several other flavours of QMD. Peilert et al. use an extension of QMD with 
additional implementation of a so-called Pauli potential PEfl . These models use a strong 
repulsive potential which depends on the distance of the particles in phase space. It is 
effective in momentum- and in configuration space and prevents two identical particles from 
coming too close in phase-space. Its parameters have been adjusted to the temperature- and 
density-dependence of the energy per particle of an ideal Fermi-gas [[^[TCj. With aid of such 
a potential selfconsistent nuclear ground states for nuclei with iV neutrons and Z protons 
as well as for infinite nuclear matter can be constructed by searching for the minimum in 
the multi-dimensional potential-energy surface of N neutrons and Z protons. The nucleons 
carry their proper Fermi-momentum, however due to the momentum-dependence of the 
Pauli-potential, their velocities (=kinetic momenta OH /dp) vanish in the ground state. This 
can be interpreted as a first approximation to antisymmetrization in finite nuclei on the two- 
body level. However one should note that for the scattering of individual nucleons a Pauli 
potential and antisymmetrisation yield different effects. Konopka et al. generalized this 
concept by treating collisions of Gaussians instead of point particles as it is the case in all 
other QMD versions. For the sake of numerical feasibility the cross section has to be assumed 
as constant. This model has been utilized for analyses of the FOPI-data at lower energies 
||71|| . The basic differences in the observables calculated with the Pauli-QMD and with 
IQMD can be explained by the use of a isotropic 41 mb cross section in the Pauli-potential 
QMD. However, one should stress that the Pauli-potential has - due to the strong momentum 
dependent potential - a different physical input as compared to all flavours discussed and 
hence a detailed comparison is not intended in this paper. 

Further studies with QMD have been done by Jaenicke et al. who replaced potentials 
and cross sections of BQMD by those calculated from a Briicker G-matrix |[72|| . Comparison 
of experimental data with this model has been performed by the FOPI-collaboration |68| . 

Lehmann and Puri extended HQMD by including a relativistic covariant propagation 
scheme of the RQMD-type. The physical inspiration of the scheme was taken from the 
RQMD -model of Sorge ]73 ], which originally was footed on the IQMD and vastly extended 
for the description of high energy collisions in the potential and the collision parts. Similarly, 
the RQMD of Lehmann and Puri is an numerical extension of the HQMD model. The 
inclusion of a covariant treatment of initialisation, Pauli-blo eking and potentials yield at 
high energies (E/A > lGeV) some differences to the 'nonrelativistic' HQMD which are 
described in |75|J76| . The relativistic and the nonrelativistic version agree at low beam 
energies. Hence this program allows a systematic investigation from very low energies to 
very high energies. The required computing time, however, is one order of magnitude higher. 

Kaon production has been intensively studied |80| using a modified version of HQMD. 
The differences are the neglect of free pions and the parametrization of the inelastic cross 
sections where only the reaction channel NN — > iVA has been employed. Differences in the 
kaon production between this version and IQMD have been discussed in 
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IV. NUMERICAL TEST AND RESULTS 



The above discussed QMD versions allow for simulations of heavy ion collisions up to 
~ 2 GeV/nucleon. Above this energy higher resonances, which are not included in the 
models under consideration become more and more important. The model gives detailed 
information about all one-body observables, such as single particle spectra, and many-body 
observables, such as particle correlations and fragment formation, on an event-by-event 
basis. Thus the structure of these theoretical data is analogous to experimental data. The 
independent development of BQMD and IQMD including different model assumptions lead 
to different results in some cases. In this section, we therefore compare several QMD versions 
with particular attention to some standard observables. It is demonstrated that most of the 
differences are related to the different treatment of the initialisation of the two colliding 
nuclei. 



1. Rapidity distributions 

A quantity, which is crucially related to the possible formation of a thermally equilibrated 
source is the rapidity distribution of baryons. 

Fig. ^ shows the rapidity distribution dN / dy of nucleons in the reaction Au+Au, b=3fm 
at 1 GeV/nucleon incident energy. IQMD (squares, dotted line), BQMD (circles, full line) 
and VUU calculations (triangles, dashed line) using a hard equation of state without mo- 
mentum dependent interaction give quite similar results. BQMD shows a slightly broader 
distribution than IQMD and VUU. As it has been already stated in ref. f4l[ the rapidity 



distribution depends strongly on the collision term and only slightly on the used nucleonic 
potentials. From this we can conclude that the hard collisions do not lead to large differences. 

The remaining differences of about 40 units at midrapidity divides up as follows (the 
statistical error of each of the midrapidity dN/dy values amount to about 10-15 units): The 
change from large to small width enhances the value by about 20 units, whereas the different 
cross sections contribute a lowering by 30 units when switching from IQMD to BQMD. The 
BQMD initialisation lowers by about 20 units and the Yukawa potential in BQMD reduces 
by about 10 units. 

At lower bombarding energies the dynamics is no longer dominated by the hard collisions 
and the nucleon potential becomes more important. At 150 MeV/nucleon, the rapidity 
distributions exhibit some larger differences between the two QMD version used. This, 
however, is due to the inclusion of Yukawa forces in BQMD. The results look more alike 
if this term is omitted in both calculations. The dN/dy at midrapidity reaches about 11% 
smaller values in BQMD, which are decomposed as follows: 7% enhancement due to the 
interaction range, 9% reduction due to the cross sections, 5% reduction due to Yukawa, and 
a slight reduction (^ 2%) due to different initialisations. 



2. Transverse flow 

Let us now focus on an observable whose investigation is strongly motivated by its 
dependence on the nuclear equation of state P, |27| , |32| (besides its dependence on the collision 
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term and on the centrality), namely the transverse flow in plane. This variable turns out 
to be extremely sensitive to a lot of parameters as we will see. The amount of transverse 
flow created in heavy ion reactions is known as a measure of the pressure built up during 
the reaction and it thus can provide information about the underlying equation of state. 

Fig. [| compares the excitation functions of flow for VUU, IQMD and BQMD with their 
default width parameters L = 8.66 fm 2 and L = 4.33 fm 2 respectively. It is found that 
VUU and IQMD show a similar behaviour with a rise of the flow up to 1 GeV incident 
energy (which is also in good agreement with experimental data) while BQMD shows rather 
weak rise of the flow. This weak rise is in disagreement with experimental flow data. The 
reasons of the differences between BQMD and IQMD shall be briefly investigated. 

The transverse flow is not only sensitive to the repulsion of the compression zone formed 
by excited nuclear matter, but also to surface properties, such as the range of the nuclear 
interaction. This quantity may be varied within the QMD approach in two different ways: 

The range of the optional Yukawa force in QMD is an adjustable parameter, it can be used 
to stabilize the width of the nuclear surface of a given density profile. The width parameter 
L of the gaussians serves as an effective interaction range as well. It should be noted that a 
change of the interaction range also changes the density gradient in inhomogeneous systems 
(this can be demonstrated by regarding the density profile of a 'box') and therefore directly 
enters into the gradient of the potential. It was found that default BQMD calculations with 
a Yukawa potential yielded a directed transverse momentum, pf r which is about 10 MeV/c 
higher than for calculations where only the Skyrme interaction is used. This is due to the 
fact that a finite range Yukawa smears out the potential gradient more than a ^-function 
and hence reduces the force in transverse direction. 

In IQMD the inclusion of Yukawa forces does not give significant effects on the nucleonic 
flow. It should however be noted, that in IQMD the range of the Yukawa force is smaller (0.4 
fm as compared to 1.5 fm of BQMD) and that actio=reactio is respected for the two-body 
interactions. 

Both models agree in the observation that a broadening of the Gaussian width L reduces 
the flow. This also corresponds to the fact that the density gradient to the high density 
region is smeared out. 

The influence of the interaction range on the flow can be studied in fig. ^ which compares 
IQMD results of the flow for L = 8.66 fm 2 and L = 4.33 fm 2 . In IQMD the default value for 
Au+Au is chosen to be L = 8.66 fm 2 . A smaller interaction range enhances the flow value 
by about 10 MeV/c at 400 AMeV and by 20 MeV/c at 1 AGeV. A further difference caused 
by the interaction range is the density of the saturation of the potential, i.e. the density 
where the potential supports maximum stability of the initial state versus vibration modes. 
For the IQMD initialisation the maximum stability is reached for L = 8.66 fm 2 at about 
p = 0.17fm" 3 and at about p = 0.15fur 3 for L = 4.33 fm 2 . 

The differences in the flow results between BQMD and IQMD motivated to search for 
parameters which might influence the flow. There are three major differences between 
BQMD and IQMD calculations: Besides the interaction range they concern the initialisation, 
the cross-sections employed and the different values for the saturation density. For a better 
comparison we changed in the following the interaction range of IQMD to L — 4.33 fm 2 . 
As a first step we investigate the dependence on the initialisation of the nuclei, which also 
includes the role of the value of po- The BQMD and IQMD initialisations differ in three 
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aspects: a) the shape of the coordinate distribution of the particles, b) the average central 
density in the nucleus and c) the limitation of the Fermi momentum to the value obtained 
by a local density approximation or the full Fermi momentum, respectively. 

The dependence on each of these differences is studied in fig. |^ which shows results of a 
HQMD calculation with inclusion of modules from BQMD and IQMD. It demonstrates that 
the very same dynamics, i.e. same forces, same cross-sections and same equations of motion 
lead to considerably varying results depending on the initial conditions chosen. 

Here IQMD-ini (squares) denotes the default IQMD initialisation with L = 4.33/m 2 . 
using hard sphere for the centroids, no constraint to the Fermi-momentum and p irii = 
Po = 0.17/m -3 . The diamonds describe a calculation with a different initial density 
Pini = 0.15/m~ 3 . We see a reduction of the flow at highest densities. This effect is known 
from hydrodynamical studies using the Rankine-Hugoniot equations. p3| . It should be noted 
that the values obtained with L = 4.33 obtained at pi n i = 0.15/m~ 3 (i.e. at maximum sta- 
bility) are comparable with the values obtained with L = 8.66 fm 2 and pi n i = 0.17/m~ 3 i.e. 
at maximum stability for the L = 8.66/m 2 case. 

The triangles denote a calculation where additionally the Fermi momenta are constrained 
by the binding energy (similar as done in BQMD). Here we see an enhancement of the flow 
at low energies. Besides the flat shape of the excitation function this constraint on the Fermi 
momentum also causes strong fluctuations of the rms radii ||. As a consequence - as we 
will see later - this introduces a considerable systematic error of the flow values. 

BQMD-ini (circles) finally denotes the BQMD default with a Wood Saxon distribution, 
the local constraints of the Fermi momentum and the saturation density of pi n i = 0.15/mr 3 . 
Since the density profile is now smeared out even more an additional reduction of the flow 
can be found at high energies. 

The composition of all three effects causes the BQMD-initialisation to yield a very flat 
excitation curve while for the IQMD-initialisation a strong dependence of the flow on the 
incident energy is observed. 

It should be noted that similar effects as reported for the flow are also found for the 
particle production. These effects are weaker but lead to the same picture. Effects that 
simulate a weaker repulsion and thus cause a weaker flow will yield an enhanced particle 
production. 

The collision term also influences the excitation function of flow. It was found that 
the collision terms of IQMD and HQMD yielded about the same values while the BQMD 
collision term causes a decrease of about 5-10 MeV/c relativ to the IQMD collision term. 

The different p d x tr values reflect themselves in different dependences of p x on the rapidity. 
Figure |6| compares the transverse momentum in the reaction plane p x (y) for the system 
Au(lAGeV)+Au at b=3fm impact parameter. We see that IQMD gives values similar to 
VUU while BQMD yields much lower flow values close to beam and target rapidity. It 
should furthermore be noted that similar effects have been found in the analysis of the flow 
out of plane. BQMD yields a less pronounced squeeze-out as compared to VUU and IQMD. 

In conclusion it is found that the description of the flow depends strongly on the de- 
tailed description of the initial state of the nuclei as well as on the interaction range. The 
constraints on the Fermi momentum as used in BQMD lower the Fermi pressure and yield 
a nearly flat excitation function of the flow which is in contradiction to current data. 
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3. Fragment production 



Another key issue for heavy ion reactions is the simultaneous production of several in- 
termediate mass fragments, i.e. clusters which are heavier than a's but considerably lighter 
than typical fission products. This phenomenon, usually referred to as multifragmentation, 
has lead to numerous speculation that this may be the signal for the occurrence of a liquid 
vapor transition in nuclear matter |[5 . 



A lot of studies involving QMD models addressed this issue as well. QMD does not 
explicitly include a phase transition, even more so, it is a non-equilibrium transport theory, 
where equilibrium need not necessarily be established to be applicable to multifragmentation 



reactions, as it is the case in statistical models for nuclear fragmentation |86| . 

It should be noted that there exist problems in describing the properties of Fermi systems 
at low temperature f46| , |8~7jj . However, it should also be noted that the fragment distributions 
obtained with QMD are in the range of the different statistical models [|53|| . The difference 
between the distributions from QMD and these models is in the same order the differences 
between those models themselves. 

QMD predicted the emission of several fragments in a single event, qualitatively similar 
to the experimental observations and at lower energies also quantitatively |5T| , ^5| , p6| , ^2| . 

A realistic description of fragmentation processes within QMD is one of the most com- 
plicated tasks. At higher bombarding energies (E > 400 MeV/nucleon) fragment formation 
is already a rather rare process. At lower energies, where multifragmentation is a major 
reaction channel (between 50 and 200 MeV/nucleon), the reaction is slowed down, which 
requires a improved accuracy of the calculation. 

One aspect is that single nuclei at rest also start to evaporate nucleons and fragments 
after 50-100 fm/c. This effect has to be minimized which sets stringent conditions on the 
stability of single nuclei. Moreover the nuclear binding needs to be properly described. 

At 50 MeV/nucleon beam energy for a symmetric system, each nucleon carries 12.5 
MeV kinetic energy in the center of mass. Together with a binding energy of about 8 
MeV per nucleon, only 4.5 MeV/nucleon are available in the center of mass. This has severe 
implications on the required accuracy of the description of ground state nuclei. If the binding 
energy is missed by only 1 MeV/nucleon, then a 22% different total energy is used in the 
calculation. At 100 MeV/nucleon this uncertainty still amounts to 6%. 

One crucial aspect, as far as the fragmentation properties of QMD are concerned, is the 
interaction range which is directly related to the width of the gaussian wavepackets. More 
extended wave packets i.e. a long interaction range leads to a smaller number of fragments. 
These fragments are somewhat heavier than those fragments from simulations with smaller 
wavepackets. This behavior has essentially two reasons: in the case of broader gaussians, 
particles in a cluster are bound to a larger number of other nucleons inside the cluster. On 
the other hand, with a smaller width the fluctuations are enhanced and an excited nucleus 
dissolves more easily. 

For example, BQMD with a Gaussian width of L = 4.33 fm 2 gives 12.7 IMFs in Au 
(150 MeV/nucleon) + Au at b=3 fm. IQMD with more extended gaussians (L = 8.66 fm 2 ) 
yields 6.6 IMFs only. It should be noted that in the present analysis the charge has not been 
regarded (especially since BQMD has no explicit charges). Therefore we used 5 < A < 19 
for the numbers obtained above. However, if we employ the same Gaussian width for both 
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models we obtain almost the same results. This can be seen in Fig. |7| where the fragment 
mass spectra have been compared for both models using both interaction ranges. 

Parameters other than the interaction range, e.g. the value of the saturation density or 
the different treatment of Fermi momenta do not affect the intermediate mass fragment mul- 
tiplicity significantly. We also find no differences on the different cross sections or potentials. 

The range of Yukawa forces do not significantly influence the mass distribution of the 
fragments as it was found for the Au(150AMeV)+Au at b=3fm. This observable shows only 
dependence on the Gaussian width. 

In conclusion it can be stated that the interaction range shows strong significance on the 
fragment production. The smaller values of L = 4.33 fm 2 used in BQMD show much better 
agreement to existing data as L = 8.66 fm 2 used in IQMD. Furthermore BQMD shows a 
better stability against particle evaporation and better binding energies. This is due to the 
constraints on the Fermi momentum as well as to the Yukawa potentials. The vibration 
modes resulting from the Fermi momentum constraints do not show strong influences on 
the fragmentation, at least in central collisions. Although the initialisation does not have 
a strong influence on the fragmentation pattern in central collisions, an initialisation which 
combines both, a proper binding energy (as it is achieved in BQMD) and a proper density 
profile (as it is done in IQMD) is preferred. This achievement is one of the main design 



goals of a new molecular dynamics scheme of the QMD type [89 



4- Particle production 

Let us now turn to the question of particle production. In the regarded energy domain 
mainly the production of pions and subthreshold kaons is of interest. 

Concerning the description of pion production the results of BQMD are not regarded, 
since it has no free pions and the cross section parametrization was not adapted to pion 
physics. Instead we will compare HQMD and IQMD. It is found that the cross section 
parametrizations of HQMD and IQMD yield very similar results. The change of the inter- 
action range changes the pion multiplicity by only 5-10%. A strong influence can, however, 
be obtained from the initialisation procedure. As an example, rapidity distribution of pions 
are displayed in Fig. ||. The two calculations differ only in the initialisation, forces and 
cross-sections are identical in both cases. The BQMD initialisation (that with the lower 
density) yields about the same shape of the pion rapidity distribution as a calculation using 
the IQMD initialisation. The absolute number of produced pions is about 30-40% larger in 
the case of the BQMD initialisation. 

The difference cannot be explained neither by Pauli blocking nor by absorption effects. It 
is found that although the calculation with the IQMD-initialisation yields higher densities, 
the calculation with a BQMD-initialisation shows higher collision numbers. This corresponds 
to the fact that up to maximum compression calculations with a BQMD-initialisation loose 
compressional energy (and thus gain kinetic energy) while calculations with an IQMD- 
initialisation gain compressional energy and loose kinetic energy. 

The arguments used in the discussion of the flow (see Figure |5|) still hold for the particle 
production. The different parts of the initialisation which cause an increase of the flow yield 
correspondingly a decrease of the pion number. For Au(l AGeV) +Au b = 3fm the change of 
the initialisation density from 0.17 to 0.15fm -3 yields (for L = 4.33fm 2 ) an increase of about 
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20 % and the change from a hard sphere initialisation to a Woods-Saxon type initialisation 
another enhancement of about 15 %. The constraints on the Fermi momentum yield no 
visible influences on the pion number at 1 AGeV energies. 

It should furthermore be stressed that all regarded models (VUU, IQMD, HQMD) per- 
form a delta decay using a lifetime which is the inverse of the mass-dependent decay-width. 
Changes concerning this description might have strong effects on the rapidity distributions. 

The study of subthreshold kaon production is motivated by the strong dependence of the 
kaon multiplicity on the nuclear eos p] , [73| , |80| , ^2| . It has been found that a hard equation 
of state yields a stronger repulsion and lower densities of the compression region than a soft 
eos ||88|| . Therefore a hard eos shows stronger flow and a smaller kaon multiplicity. The pion 
multiplicity shows only slight dependences on the eos since both equations of state yield 
about the same compression densities PJ58|. 

Kaon production have not been studied within the BQMD model but with an upgrade 
version called QMDRKNC This version does not include free pions, therefore the 

deltas have an infinite lifetime. A first comparison of this version with IQMD has been 
presented in |f75 |. The kaon numbers obtained with both models agree for Au (1 AGeV) 
+Au within 20-30%. HQMD with L = 4.33fm 2 and BQMD initialisation yields about the 
same values. HQMD differs mainly from QMDRKNC in the lifetime of the deltas. The 
similar multiplicities of IQMD and QMDRKNC are a result of counterbalancing effects 
which will be briefly discussed: 

The infinite lifetime of the delta in BQMD and QMDRKNC causes an enhancement of 
about 10 - 20% of the kaon number when compared to HQMD. This is due to the dominance 
of the channel nucleon+delta — > nucleon+hyperon +kaon for the subthreshold production. 
An infinite lifetime enhances the possibilities for nucleon-delta collisions. Similar numbers 
have been found when comparing default IQMD with an calculation with infinite delta 
lifetimes. 

The BQMD initialisation yields an enhancement of the kaon production by about 20 - 
30 % as compared to the IQMD initialisation. The reason is presumingly similar to that for 
the pion production. The BQMD initialisation allows higher kinetic energies of the nucleons 
in the compressed state. 

The choice of a short Gaussian width (L = 4.33fm 2 , BQMD default) causes a reduction 
of the kaon number of about 30 % as compared to a calculation with L = 8.66fm 2 which is 
the IQMD default for Au. The reason for this may be connected to the argument used for 
explaining the enhancement of flow when using a short width. The density gradient gets 
steeper when the interaction range decreases. This simulates a stronger repulsion of the 
compressed nuclear matter. 



5. Initialisation and Stability 

One of the seminal problems of the simulation of a heavy ion reaction is the proper 
description of ground state nuclei. One cannot expect that a reaction is reproduced properly 
if projectile and target do not have the observed properties, in particular the proper ground 
state density. 

As we have stressed several times before the choice of the initial condition is crucial 
for a proper description of various phenomena. Fermi-momenta treated semiclassically as a 
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random motion of nucleons inside a nucleus induces significant fluctuations of the density 
profile, if the motion of a single nucleus is followed for some time. 

Fig. |9] shows a time evolution of the root mean square radii of single Au nuclei in 
coordinate and momentum space. Concerning the coordinate space we observe an expansion 
mode in IQMD and an oscillation mode in BQMD. IQMD shows best stability if potentials 
corresponding to a hard eos is used whereas BQMD shows best stability if a soft eos is used. 
The rms radii obtained with BQMD soft eos correspond to the results presented in M . The 
other eos yield larger fluctuations for both programs. This stresses once more the fact that 
a semiclassical approach can be optimized to a desired purpose but on the cost of other 
applications. 

The fluctuations of the rms radii in momentum space demonstrate that the potentials are 
not saturated in the given initialisation. The system converts potential energy into kinetic 
energy and vice versa to equilibrate the system. This conversion of energy was already 
addressed in the previous subsection (particle production). The initialisations of IQMD and 
BQMD yield different pion numbers due to different balances of kinetic and potential energy. 

The counterbalancing parts are the kinetic pressure which causes an expansion and the 
density dependence of the potentials which may cause attraction for low densities and re- 
pulsion for high densities. 

For this it may be interesting to regard the mean density of the system which is the 
mean value of the density of each particle averaged over all particles. It should be noted 
that this value may be sensible to density fluctuations in the center which do only slightly 
effect the root mean square radius. 

A time evolution of single nuclei yields for the Au case (hard eos) mean densities (density 
per particle averaged over all particles) changing between about p = 0.14 and 0.15 fm -3 for 
IQMD and between p = 0.11 and 0.18 fm -3 for BQMD. For smaller systems (on a time 
period of about 60 fm/c) the stability of IQMD gets smaller, to e.g. a range of p = 0.14 to 
0.17 fm -3 for a Nb nucleus in IQMD, p = 0.16 to 0.19 fm~ 3 for a Ca nucleus in IQMD and 
p = 0.16 to 0.21 fm -3 for a Ne nucleus in IQMD, while in BQMD the fluctuations remain 
constant (p = 0.1 — 0.17fm~ 3 ) for all three regarded systems. It should also be noted that 
these fluctuations increase in IQMD if a soft eos is used (about p = 0.11 — 0.15 for the Au 
case) and decrease in BQMD (to about the same values). 

Let us now examine the density profile of a single Au nucleus. Fig. [10] displays the time 
evolution of the density profile within BQMD. We observe a change in the center as well as 
at the surface. In a considerably large volume around the center (r < 5 fm) the change of 
density with time induces changes of the compressional energy in a heavy ion reaction. The 
weakening of the surface causes an increase of the rms radius and it therefore modifies the 
total interaction cross section as well as the probability of e.g. particle production processes 
in particular in peripheral collisions. Without Yukawa forces these fluctuations are even 
larger. The reason for these fluctuations are the lack of pressure built up by the Fermi 
momentum when the nucleus gets compressed. This can be verified by initializing HQMD 
with the full Fermi momentum but otherwise as above. 

The time evolution of the density profile in IQMD is displayed in Fig.|TT] for a gold nucleus 
with the default width of L = 8.66 fm 2 . We find strong fluctuations at r = but a stable 
shape at the surface and a rather stable rms radius. We also find that the shape of the 
r 2 p{r) distribution shows better stability for the L = 8.66 fm 2 case than for L = 4.33 fm 2 . 
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This also motivated the choice of L in IQMD. 

It has now to be tested whether these fluctuations of the density profile cause uncer- 
tainties in the determination of observables. This is tested by changing the initialisation 
distance d (with respect to the minimum distance for a head-on collision) which is defined 
as (R are the radii of the nuclei): 

d = ^(center of proj.) — ^(center of target) — i?(proj.) — i?(targ.) 

By changing this distance we allow the nuclei to change their profile according to the internal 
forces before they come into nuclear contact. For the ideal case this technical parameter 
should have no influence on the observables. In reality the observables depend on it, however, 
in most cases weakly. 

A strong influence of the density fluctuations on observables has been found particularly 
within BQMD for the collective sideward flow. Depending on the initially chosen distance the 
total directed transverse momentum transfer in the reaction Au (1 GeV/nucleon, b=3fm) 
+ Au varies strongly as it can be seen in fig|T2|. Although the absolute magnitude of 
the flow depends on whether a Yukawa interaction is employed or not, this variance is 
observed in both cases. IQMD, however, run with the default parameters shows a much 
weaker dependence on this technical parameter. In the calculation the flow varies as a 
function on the initialisation distance only by about 10%. An IQMD calculation with 
L = 4.33 fm 2 shows less stability. The flow values are decreasing with initialisation distance. 
This corresponds to the effect that the density of maximum stability (p = 0.15 fm -3 for 
L = 4.33 fm 2 ) is not equal to the initialisation density pj„j = 0.17fm -3 . This was originally 
the motivation for the use of a system dependent width in IQMD. A BQMD calculation with 
L = 8.66, however, still shows strong fluctuations. The fluctuations in BQMD decrease if 
the Fermi momentum of the initialized nuclei is increased. However, in this case the binding 
energy of the nucleus becomes smaller and spurious particle evaporation may be effected. 

For the observables for Au(lGeV)+Au discussed in this paper we find the following max- 
imum deviation (in a range of initialisation distances between and d max = 13 fm) from 
the default values: for the IQMD initialisation 10% concerning flow and kaon multiplic- 
ity and 8% concerning pion multiplicity and for the BQMD initialisation 70% concerning 
flow, 45% concerning kaon and 35% concerning pion multiplicity. For the fragmentation of 
Au(150Mev)+Au b=3fm both models yield (in their default modes) about 8-15% deviation 
in the number of IMFs ( a rise from 6.6 to about 7.1 for IQMD and a fall from 12.7 to about 
10.8 for BQMD). 

It should be noted that the reported errors also include about 5-10% statistical fluctua- 
tions and that for initialisation distances larger than 13 fm the deviations may still increase 
for some calculations. The value of 13 fm has been chosen since it is the difference in the 
effective distance to the first reaction point between a central and a very peripheral collisions. 

In any case, these fluctuations cause an additional systematic error which has to be added 
to the statistical one if one compares with data. Because these fluctuations are stronger using 
a BQMD initialisation the systematical error is larger there. 

In conclusion we find that both models show fluctuations of the density profile at r = 0. 
However IQMD shows in its default mode a rather stable r 2 p(r) shape while for BQMD the 
maximum of r 2 p(r) changes in time. This yield artificial vibration modes which influence 
the stability of the nucleus and therefore cause systematic errors which in most cases are 
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stronger in BQMD than in IQMD. Especially the discussion of dynamical variables like flow 
and particle production within BQMD has to be regarded very cautiously. 

V. SUMMARY AND CONCLUDING REMARKS 

We have compared different realizations of the Quantum Molecular Dynamics model. 
The different realizations differ in certain input variables as they are comparatively presented 
in table |IV} Some of these parameters are purely technical, some are physical. The latter 
parameters are constrained by experimental observations but are not completely fixed. If the 
same parameters are employed the result of the different programs are identical in between 
the error. This is a remarkable achievement in view of the several thousand program lines 
of each of these programs. 

The HQMD realization allows for the first time to compare in detail the influence of 
different inputs on the different observables. The most important input is the choice of the 
Gaussian width L or in other words the interaction range of the nuclear potential. A change 
of the interaction range causes differences in the density profile of a ground state nucleus and 
in the strength of the density gradient. Thus a smaller interaction range yields an enhanced 
flow, enhanced fragment production and (which was not shown) a reduced numbers of pions 
(slight changes) and kaons (larger changes). The interaction range determines the surface 
properties of the nuclei as well as their binding energy. Only in infinite nuclear matter the 
binding energy is independent of this quantity. The experimental value of these observables 
allow to fix the range of possible values. 

The choice of the cross section employed in HQMD, BQMD or IQMD yield slight changes 
in the flow (and also slight differences in the stopping) but has no influences on the frag- 
mentation. Pion production in HQMD and IQMD are comparable. 

An important factor is the choice of the ground state description. This choice effects the 
results in flow and particle production, but does not influence the fragmentation pattern. 
The flow at high energies is found to be stronger if the initialisation of the nucleus is more 
compact. At low energies these differences vanish, instead the average Fermi momentum 
becomes important. The Fermi momentum corresponding to (infinite) nuclear matter are 
necessary to stabilize the nucleus against artificial vibrations and yield better agreement 
with experimental flow data. At the same time a large Fermi momentum lowers the bind- 
ing energy. The binding energy of heavy nuclei is reproduced if local Fermi momenta are 
employed. 

As a conclusion we find that there are variables which are very robust against a change 
of the technical or physical parameters, e.g. the rapidity distribution. Others like the frag- 
mentation pattern depend on the range of the interaction only. Other observables, like the 
directed flow have a very strong dependence on many details of the calculation and slight 
differences between the different QMD flavours yield large differences in this observable. 

In the BQMD proper binding energies have been achieved on the cost of large dy- 
namical fluctuations in the initial state and by a moderate energy non- conservation (« 
2 MeV/nucleon for single nuclei) on an event by event basis. The IQMD-model offers rather 
stable density distributions and good energy conservation, however for the price of nucleon 
evaporation and and improper binding energies {E^nd ~ 4 — 5 MeV/nucleon for heavy nuclei 
instead of 8 MeV/nucleon). 
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The choice between a parametrization which yields the proper ground state energy of 
projectile and target and that which yield the necessary Fermi momentum to obtain the 
observed flow is not satisfying. Therefore work is in progress to modify the bare interaction 
between the nucleons inside the nucleus in a way which allows to obtain both at the same 
time. In particular one has to account for the peculiar dependencies between the various 
parts of the model. Initialisation, propagation, hard collisions, and Pauli-blocking cannot be 
treated independently from each other. Forces and cross-sections are connected. The specific 
choice of the saturation density, and thus also that of the central density of heavy nuclei 
influences the physical output. These aspects are part of the effort to obtain a new unified 
QMD scheme which covers the energy range between 25 MeV/nucleon and 200 GeV/nucleon. 
The aim of this new model will be to cover the best possibilities for a reliable description of 
the different aspects of heavy ion collisions fl89|. 
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TABLES 



a (MeV) 



P (MeV) 



7 



5 (MeV) 



feet) 



S 

SM 
H 

HM 



-356 
-390 
-124 
-130 



303 
320 
71 
59 



1.17 
1.14 
2.00 
2.09 



1.57 
1.57 



500 
500 



TABLE I. Parameter sets for the nuclear equation of state used in the QMD model. S and 
H refer to the soft and hard equations of state, M refers to the inclusion of momentum dependent 
interaction. 



x = y/s (GeV) a (fm) b 



2.104 - 


- 2.12 


294.6 (x - 2.014) 2 - 578 


19.71 (x - 2.014) 1 - 551 


2.12 - 


2.43 


0.01224 
(x-2.225) 2 +0.004112 


19.71 (x - 2.014) 1 - 551 


2.43 - 


4.50 


(2.343/x) 43 - 17 


33.41 arctan(0.5404 (x - 2.146) - 9784 ) 



TABLE II. a(s) and b(s) as functions of the cm. energy. 



x = \Ap(CM)\/l GeV o e i (mb) for pp, nn a e \ (mb) for pn 



x < 0.8 


23.5 + 1000- (0.7 -x) 4 


33 + 196- |0.95 -x| 2 - 5 


0.8 < x < 2 


1250/(x + 50) - 4-(x-1.3) 2 


31/v^ 


2 < x 


77/(x + 1.5) 


77/(x + 1.5) 



TABLE III. Elastic cross section parametrization used in QMDRKNC and HQMD as a func- 
tion of the relative momentum in the CM-frame 



Input 


BQMD 


IQMD 


HQMD 


Initialisation 


Wood- Saxon 


hard sphere 


both 


Init. distance 


3fm 


Ofm 


3fm 


Gaussian width 


4.33 fm 2 


8.66 fm 2 


4.33 fm 2 


Coulomb forces 


Z p = Z n 


Zp = 1, Z n = 


Zp = Z n 


Yukawa forces 


L=1.5 fm 


L=0.4 fm 


none 


Yukawa adjust. 


ti(i) = t 1 -E j u? j uk /p 


h = 

h-n* Vf uk * L Yuk /L 3 / 2 


none 


actio = reactio 


on the average 


'exact' 


'exact' 


asymmetrie forces 


none 




none 


forces on ir 


no 7T 


Coulomb 


no force 


cross sections 


Cugnon 1981 


VUU 1986 


Cugnon 1989 


particles 


N, A 


N, A, 7T 


N, A, N*,ir 



TABLE IV. Comparison of the different realizations BQMD, HQMD and IQMD concerning 



the different ingredients of the inputs. 
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FIGURES 



FIG. 1. The elastic and inelastic cross sections for proton-proton (pp) and proton- neutron 
(pn) used in IQMD. The neutron-neutron cross section is assumed to equal to the pp case. The 
total cross section is equal to the sum of elastic and inelastic cross section. 

FIG. 2. Rapidity distributions dN/dy of nucleons in the reaction Au(l GeV/nucleon)+Au, 
b=3fm for VUU, IQMD and BQMD. In all calculations a hard equation of state without momentum 
dependent interactions has been used. 

FIG. 3. Excitation function of the system Au+Au at b = 3fm impact parameter obtained 
with BQMD, IQMD and VUU in their default versions. 

FIG. 4. Excitation function of the system Au+Au at b = 3fm impact parameter obtained 
with IQMD using the width of L = 4.33 fm 2 and L = 8.66 fm 2 . 

FIG. 5. Excitation function of the system Au+Au at b = 3fm impact parameter obtained 
with HQMD (default collision term) using the BQMD initialisation, a hard sphere initialisation 
with reduced Fermi momentum and with full Fermi momentum and the IQMD initialisation. 

FIG. 6. Comparison of the transverse flow p x (y) of nucleons in the reaction Au(lAGeV)+Au 
b=3fm for VUU, IQMD and BQMD, all using a hard equation of state without momentum depen- 
dent interactions. 

FIG. 7. Fragment mass distribution obtained by BQMD and IQMD for the system 
Au(150MeV) +Au b=3fm, both with L = 4.33 fm 2 and L = 8.66 fm 2 . 

FIG. 8. Rapidity distributions of pions for Au(l GeV/nucleon)+Au reactions at b=3fm ob- 
tained with QMD employing the BQMD initialisation (circles) and the IQMD initialisation (dia- 
monds). 

FIG. 9. Time evolution of the root mean square radii of a single Au nucleus in coordinate and 
momentum space obtained with IQMD and BQMD using a hard and a soft eos. 

FIG. 10. Time evolution of the density profiles p{r) obtained for a Au nucleus initialized with 
BQMD (with Yukawa) using the width of L = 4.33 fm 2 . 

FIG. 11. Time evolution of the density profiles p(r) obtained for a Au nucleus initialized with 
IQMD using the width of L = 8.66 fm 2 . 

FIG. 12. Dependence of the flow obtained from the system Au(l AGeV)+Au at b=3fm on 
the initialisation distance using IQMD with L = 8.66 fm 2 and BQMD with L = 4.33 fm 2 . 
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Figure 
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Au+Au b=3 fm, hard eos 
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Au+Au b=3 fm, hard eos 
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Au(lAGeV)+Au b=3fm 
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Au(150AMeV)+Au b=3 fm, hard eos 
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single Au nucleus 
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Au nucleus in BQMD (L=4.33) 
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Au nucleus in IQMD (L=8.66) 
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Figure J5 
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Au(l AGeV)+Au b=3 fm, hard eos 
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Figure J2 



39 



